from coopr.pyomo import *

#
# Model
#
model = AbstractModel()

model.I = Set()
model.J = Set()

model.A = Param(model.J, model.I)
model.b = Param(model.J)
model.c = Param(model.I)
model.p = Param(model.J) #probability for scenarios

model.alpha = Param()    #required chance
model.bigm  = Param(model.J)

model.X = Var(model.I, within=PositiveReals)
model.Y = Var(model.J, within=Binary)

def knapsack_rule(mod):
    return sum(mod.p[j] * mod.Y[j] for j in mod.J) >= mod.alpha
model.KnapsackRule = Constraint(rule=knapsack_rule)

def chance_rule(mod, j):
    lhs = sum(mod.A[j, i] * mod.X[i] for i in mod.I) - mod.b[j]
    rhs = mod.bigm[j] * (mod.Y[j] - 1)
    return lhs >= rhs
model.ChanceRule = Constraint(model.J, rule=chance_rule)

def total_cost_rule(mod):
    return sum(mod.X[i] * mod.c[i] for i in mod.I)
model.TotalCost = Objective(rule = total_cost_rule, sense=minimize)


